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"IHEIUduc L ion ;  )  In  Lli  and  Rosenblatt  (1982)  a  deconvolution  scheme  for 


nonGaussian  linear  processes  making  use  of  third  order  moments  (or 
spectra)  was  presented.  This  is  appropriate  for  such  processes  with 
nonzero  third  order  central  moments.  However,  if  the  third  order 
momenta  are  zero  (this  could  happen  in  the  case  of  symmetric  distribu¬ 
tions)  it  is  appropriate  to  look  for  a  fourth  order  technique  that 


would  be  effective.  Such  a  scheme  is  presented/  and  discussed  in  this 


paper  together  with  some  illustrative  examples! 


We  give  a  brief  sketch  of  the  theoretical  background.  Let  vfc, 
t  -  ...,-1,0,1,...  be  Independent,  identically  distributed  random 
variables  with  mean  zero  and  variance  one.  Consider  a  sequence  of 


real  constants  {a^}  with 
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Let  {xfc}  be  the  linear  process 


(1) 


Yt-j 


Introduce  the  z-transform  a(z)  ■  \  a.z^  corresponding  to  the  process 

J  3 

{xt>.  He  should  like  to  estimate  a(e~^)  (X  real)  from  observations 
only  on  the  process  {x^}  and  use  this  estimate  to  deconvolve  the  pro¬ 
cess  {xfc}  and  estimate  the  £t's.  As  noted  in  Rosenblatt  (1980),  in 
the  Gaussian  case  one  can  only  estimate  the  modulus  of  a(e  **),  and 
it  is  only  in  the  nonGausslan  case  that  one  can  also  estimate  the  argu¬ 
ment  of  o(e  **).  Of  course,  the  spectral  density  of  {xfc}  is 


f(X)  -  |a (e-1*))2  . 

In  some  geophysical  contexts  a  nonGausslan  model  like  that  discussed 
here  has  been  proposed.  A  basic  concern  is  that  of  dec  evolution,  esti¬ 
mating  the  a^'s  and  v^'s.  A  discussion  of  such  questions  with  some  of 
the  geophysical  background  can  be  found  in  Donoho  (1981),  Godfrey  and 
Rocca  (1981) ,  and  Wiggins  (1978) . 

The  following  lemma  was  proved  in  Lii  and  Rosenblatt  (1982)  with  the 
type  of  argument  suggested  in  Rosenblatt  ( 1980) . 


Lemma.  Consider  a  nonGauasian  linear  process  (xfc)  (see  (1))  with 
the  independent  random  variables  having  all  their  moments  finite .  Let 


^|J!  !«jf  <  " 

and  assume  a(e_iX)  *  0  for  all  X.  The  function  a(e"iX)  can  then 
be  identified  in  terms  of  observations  on  only  {xf}  up  to  an 
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indeterminate  integer  a  in  a  factor  and  an  indeterminate  sign 
of  a(l)  “  \  a^.  For  this  result  it  ie  sufficient  to  have  some  finite 
moment  of  order  k  >  2  with  cumulant  +  0. 

The  ktb  order  cumulant  spectral  density  of  the  process  {x^}  can 
be  seen  to  be 


-iX, 


-iX. 


b  (X......X.  .)  - t— r  a(e  *)...  (e  K"A)o(e 

k  1  k  1  <2ir)k_1 


i(X,  +•••  +X.  ,) 


)  . 


If  one  sets 


htt)  -  .rg  {aC1*) 


then  it  can  be  shown  that 


h(Xt)  +•• 


•  +h(Xk_1)  -  h(X1  +  -»» 


“  ars  [{  1«S}|}  y*  VX1 . Xk-1)] 

The  case  In  which  k  *  4  is  of  particular  interest  to  us.  From  this 
point  on  we  will  often  delete  the  subscript  k  »  4  but  understand  that 
we  are  dealing  with  4tb  order  cumulant  spectral  estimates.  The  rela¬ 
tion 
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■*(u)  -  h'*(0))du  +  cX  *  h^(X)  +  cX  . 
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Since  h(ir)  must  be  an  integral  multiple  of  w  (because  the  c^'s 
real)  we  can  rewrite  ( 2 )  as 

h.  (ir) 

h*  (X)  -  h,  (X)  -  — - X  +  aX 

l  ir 


with  a  an  indeterminate  Integer.  Further 


h'(0)  -  h" (X)  -  lim  jr  (h(X)  +  2h(A)  -  h(X+2A)> 
A+0  L 


and 


h(X)  +  2h(A)  -  h(X+2A)  -  arg  (b(X,A,A)} 

up  to  a  sign.  We  shall  consider  the  question  of  estimating  hx(X) 
Set  A  -  A(n) ,  (2k+l)A  -  X  and  consider  A  ■  A(n)  0  as  n  •. 
Clearly  b (0,0,0)  is  positive.  Notice  that 


hx(X) 


-  h(X)  -  h'(0)X 

«h((2k+l)A)  (2k+l)A 

k-1 


-2  (h((2J+l)A)  +  2h(A)  -  h((2j+3)A)> 
J-0 


k-1 

E 


arg  b((2J+l)A,A,A) 


j-0 
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are 


This  suggests  taking 


(3) 


H  (X)  -  arg  b((2j+l)A,A,A) 

j-0 

as  a  possible  estimate  of  h^(X).  We  shall  assume  that  Qb  is  a 
consistent  sequence  of  estimates  of  the  fourth  order  cumulant  spectral 
density.  Conditions  for  existence  of  such  a  sequence  of  estimators  can 
be  found  In  Brillinger  and  Rosenblatt  (1967).  Here  n  denotes  the  sample 
size.  Take 


6  (X.y.n)  ■  arctan  (Im  b(X,ytn)/Re  b(X,y,n)) 
u  n  n 

as  an  estimate  of 


8(X,y,n)  •  erg  b(X,y,n)  . 

Then  just  as  in  Lli  and  Rosenblatt  (1982)  one  has 

en(x,u,n)  -  e(x,y,n) 

-  -  lg  b(-^4  (Re  b(X,y,n)  -  Re  b(X,y,n)> 

|b(Xty ,n) \  ° 

+  {lB  b(x,ytn)  -  Im  b(X,y,n)} 

|b(X,y,n)r 

+  <^(nb(X,y,n)  -  b(X,y,n))  . 

-IX  2 

Suppose  Hr(X)  Is  taken  as  an  estimate  of  h^(X).  If  a(e  )  €  C  y 

the  weight  function  of  nb  Is  symmetric  and  bandlimited  with  bandwidth 
3 

A,  A(n)  +  0,  A  a  +  ■  es  n  •,  then 
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«n(X)  -  h^A)  -  Rn(A)  +  op(Hn(A)  -  hj (A ) ) 

where 

R  (A)  -Y  ria  b((2J-H)AiAiAl  {Re  b(2j+l)A,A,A)  -  Re  b((2j+l)A,A,A)} 
n  llb((2j+l)AfA,A)|2 

(4)  -  Re.  b<<2J+il)AiA.iA).  {In  b«2j+l)A,A,A)  -  la  b((2j+l)A,A,A)}l 

|b((2j+l)A,A,A)r  J 

One  can  show  that 


E  R  (6) 
n 


—jf 


b(ul’u2,U3)|u2-u3-0  du  4 


+  o(A) 


where  the  A^k  are  the  moments 


jk 


■  /"kUk  W(ul,«2.u3)d»1 


du2du3 


and  Du  is  the  partial  derivative  with  respect  to  .  Further 


2  J»in(A,p) 


cov 


(R(A),Rn(|i))^^  /  {f2(u)/|b(u,0,0)|2}du 

£  An  •'0 

/WVv, 

f  2 

[»  (A,p)  I  W  (u,v,w)dudvdw 


,w)dudvdw  . 


4* 
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Here  W  la  the  standardized  weight  function  of  the  fourth  order  cumulant 
spectral  density  estimate. 
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Computational  methods.  We  consider  computational  schemes  for  computing 

Hn(*)  1“  (  3)  as  an  estimate  of  h^X).  Given  a  sample  (xt)  of  size 

n  *  mN  break  up  the  sample  into  m  disjoint  subsections  of  equal 

length  N  so  that  the  variance  of  the  trispectral  estimate  (estimate 

of  the  cumulant  spectral  density  of  fourth  order)  of  each  section  is 

not  too  large.  This  is  particularly  important  in  the  trispectral  case 

since  the  variance  of  the  fourth  order  periodogram  used  in  constructing 

2 

trispectral  estimates  is  proportional  to  N  .  Then  choose  a  grid  of 
points  Xj  “  (2j+l)A  in  (0,2ir),  j  ■  0,l,...tM,  A*  2irL/N  for  a  suitable 
integer  L.  Form  the  trispectral  estimates  b(X  ,A,A),  s  ■  l,...,m, 

®  j 

for  each  of  the  m  sections  of  length  N  by  using  a  weight  function  of 
bandwidth  2A,A,A  in  each  component  since  X ^  -  A  -  2A.  Average  the 
estimates  obtained  from  each  of  the  subsections  of  length  N  at  X^  to 

A 

obtain  the  final  estimate  nb(Xj,A,A).  Compute  0n(Aj)  “ 

arg  {  b(X  ,A,A)}  and  form 
U  J 

i-1 

H  (A  )  ■  -2e„ <*4>»  1  "  1*2, •  ,M+1  . 

3  j-0  "  3 

We  set  H  (0)  ■  0  since  h(0)  ■  0  and  estimate  H  (X_)  -  H  (A)  by 
n  no  n 

an  interpolation  between  0  and  H  (X.)  ■  H  (3A) .  Then  coefficient 

n  x  n 

in  the  trigonometric  expansion  of  a(e~^)  can  be  estimated  by 


(5) 


■H 


dX 


(2ir  f^))* 


HU) 

exp  U(H  (X.)  -  — - X.)  +  ikX.} 

n  J  if  j  j 


where  fn(X^)  is  a  consistent  estimate  of  the  spectral  density  f(X) 

of  {x  }.  The  spectral  density  estimate  f  (X  )  is  formed  as  follows: 
^  “  J 

Form  smoothed  periodograms  with  bandwidth  A^  5  2A  from  each  of  the  m 

subsections  of  length  N  and  average  the  m  smoothed  periodograms  to 

get  f  (X.).  f  (o)  is  estimated  by  extrapolation.  Presumably  one 
n  j  ** 

could  Improve  formula  (  5)  by  using  a  more  refined  approximation  to 
the  Integral  based  on  the  trapezoidal  rule  or  Simpson's  rule.  Also 
an  extrapolation  procedure  could  be  used  at  the  end  points  since 
0,A,3A, . . .  are  not  equally  spaced. 


We  now  describe  an  alterative  procedure  for  estimating  h^(X) . 
Note  that 


k-1 

2J  arg  b(jA,A,A) 

J-l 


k-1 

■r 


(h(JA)  +  2h(A)  -  h(jA+2A)> 


j-l 

2[kh(A)  -  h(kA) ]  +  B 


with 


B  -  h(2A)  -  h(A)  +  h(kA)  -  h((k+l)A)  . 


Thus  If  X  -  kA 


(6) 


h^A)  -  h(X)  -  h-(0)X 
k-1 


a  -  4  2  arg  b(jA,A,A)  -  ~  B 

j-1  2 


If  A  Is  small  we  would  expect  B  to  be  small  also.  This  suggests  that 
a  plausible  estimate  of  h^(X)  could  be  given  by 


Gn(X) 


k-1 

-  j  5-)  arg  b(jA,A,A) 
j-1  n 


The  estimate  Gn(A)  may  have  an  additional  bias  relative  to  the  estimate 

Hq(A)  because  of  the  term  -  yB  in  (  6  ) .  However,  a  full  comparison 

of  the  two  estimates  is  difficult  to  make.  There  are  advantages  and 

disadvantages  to  each.  The  estimate  actually  used  in  the  computational 

illustrations  discussed  later  is  G  (X). 

n 

To  deconvolve  the  observed  signal  {x£>  and  obtain  estimates  of 
(vt),  we  form 

vt  -  a  A(L)xfc 

where  L  is  the  backward  shift  operator.  When  o(e”^)  is  one  sided 
polynomial  of  order  q  (this  corresponds  to  {xfc}  a  moving  average  of 
order  q)  methods  using  a  partial  fraction  expansion  of  a  *(L)  by  com¬ 
puting  the  roots  of  o(z)  are  described  in  Lil  and  Rosenblatt  0.982  ) . 


To  avoid  finding  an  appropriate  finite  parameter  model  for  {xfc}  and 
dealing  with  the  sensitivity  of  root  location  in  terms  of  their  depend¬ 
ence  on  coefficients,  we  note  that  one  can  find  the  deconvolution  weights 
by  inverting  a(e-1A)  directly.  Let  b(e-1X)  -  o(e_iX) .  Then  the  coef¬ 
ficient  bfc  in  the  expansion 


sf  **)  -23  \  e“d 


can  be  computed  by  using 


X2ir 

-4  GOO 

[2*  fn(X)]  *  exp  (l(-Gn(X)  +-\-\  +  kX)} 

J 

(2'  W —  <«*«■<»!>  +  ^  +  *V’ 

k  “  ...,-1,0,1,...  .  Usually  we  find  uitable  integers  and  k^  and 

use  the  real  part  of  b^  for  k  -  k^, ... ,k2  as  deconvolution  weights 
since  we  are  dealing  with  a  real  process.  In  the  examples  discussed 


below  the  choice  was 


k^  *  -  9  and  k2  -  9, 


Examples.  A  few  simple  examples  are  presented  here  to  illustrate  the 
computational  procedures.  The  model  considered  is 


t  -  1 


vt  -  (v£  -  v')/s 


I 


I 


I 

I 


640 


V*  ■  /  ^  V* 


/640 


2 

8 


£ 

t-i 


(v'  -  v')2/640 


and  the  v^'s  are  independent  and  identically  distributed.  The  general 
computational  set  up  is  the  same  as  that  in  Lii  and  Rosenblatt  (1982). 
All  the  examples  deal  with  schemes  generated  with  coefficients  (and 
roots)  as  specified  in  Table  1. 


Table  1.  Coefficients  and  roots  for  four  cases 


Coefficients 

Roots 

Case 

ao 

°1 

°2 

rl 

r2 

1 

1.0 

-  .833 

0.167 

2.0 

3.0 

2 

1.0 

-2.333 

0.667 

0.5 

3.0 

3 

1.0 

-3.50 

1.50 

2.0 

.333 

4 

1.0 

-5.0 

5 

6.0 

0.5 

.333 

In  the  first  set  of  examples  (four)  v'  is  the  exponential  distribution 
with  parameter  1  generated  by  GGEXN  in  IMSL.  Although  the  third  order 
cumulant  of  v'  is  nonzero,  the  fourth  order  technique  considered  in 
this  paper  can  be  U8ed.  Table  2  compares  the  estimated  coefficients 
in  each  of  the  cases  as  computed  by  third  and  fourth  order  techniques. 


Table  2 


Third  order 


a. 


-0.661 

-2.132 

-3.32 

-3.23 


* 

°2  ! 

A 

°0 

.05  i 
.796  l 
1.17 
6.56 

.6955 

1.011 

1.44 

.805 

Fourth  order 


a. 


-  .747 
-2.043 
-3.25 
-4.456 


The  deconvolution  of  case  2  using  the  third  order  method  is  shown  in 
Figure  la.  This  can  be  compared  with  deconvolution  by  the  fourth  order 
method  which  is  given  in  Figure  lb.  Both  deconvolutions  in  Figure  1 

involved  computation  of  roots.  The  mean  square  errors  of  v  -  vfc 
for  the  third  order  and  fourth  order  methods  were  .045  and  .094  re¬ 
spectively. 

In  the  second  set  of  examples,  the  v'  distribution  was  the  sym¬ 
metric  double  exponential.  The  estimated  coefficients  for  the  four 
cases,  using  a  fourth  order  method,  are  given  in  Table  3. 

Table  3 


mm 

A 

°1 

A 

°2 

-  .3886 

-  .1 

04 

.8835 

-2.121 

.8 

02 

1.874 

-2.544 

1.1 

53 

1.805 

-3.22 

3.8 

65 

The  deconvolution  for  case  2  using  location  of  roots  is  given  in  Figure  2 
The  mean  square  error  of  v^  -  v^  is  .01559  while  the  variance  of  vfc  is 
25. 
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In  the  last  set  of  examples  v '  has  a  symmetric  Pareto  distribution. 
First  uniform  random  numbers  (on  the  interval  (0,1))  are  generated 

by  GGUW  in  IMSL.  Then  the  transformation  y^  ■  (U^)”*^  is  used  to 
obtain  random  numbers  having  a  Pareto  distribution  with  density  f(y)  = 

5y  y  >  1.  The  v^'s  are  obtained  by  randomly  changing  the  sign  of 
y^  -  1  with  probability  .5.  The  estimated  coefficients  in  the  four 
cases  using  a  fourth  order  procedure  are  given  in  Table  4. 


Table  4 


Case 

a0 

&1 

*2 

1 

1.040 

-  .6566 

.1919 

2 

1.037 

-2.095 

.5926 

3 

1.265 

-3.32 

1.252 

4 

1.599 

-5.158 

4.799 

Figure  3m  gives  the  result  of  deconvolution  of  case  2  using  computation 
of  roots.  Figure  3b  gives  the  result  of  direct  deconvolution  in  case 
2.  The  mean  square  error  of  vt  -  vt  is  .0059  and  .0011  for  the  first 
and  second  deconvolution  procedures  respectively. 

Comments  on  computation.  A  decision  as  to  when  to  use  a  third  or  fourth 
order  deconvolution  procedure  could  be  based  on  estimates  of  third  and 
fourth  order  cumulants.  A  larger  estimate  (in  absolute  magnitude)  for 
a  specific  cumulant  would  suggest  that  one  could  with  some  condifence 
prefer  using  the  deconvolution  procedure  of  the  same  order.  Of  course, 
if  the  cumulants  were  too  small  in  magnitude  there  wouldn't  be  much 
point  in  attempting  the  deconvolution. 
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The  sample  size  used  in  the  illustrative  computations  is  640.  In  the 
ordinary  usage  this  would  be  thought  of  as  a  large  sample.  One  thing  that 
is  apparent  is  the  relative  effectiveness  of  the  deconvolution  procedure 
Independent  of  the  tail  behavior  of  the  vfc  distribution.  But  one  can 
say  more.  In  a  certain  sense  the  sample  size  640  is  moderate  (perhaps 
even  small).  Suppose  we  look  at  the  question  of  estimating  the  third 
and  fourth  central  moments  when  one  has  a  sample  size  of  n  observations. 
The  first  order  expressions  for  the  variances  of  the  standard  estimates 
of  third  and  fourth  central  moments  are 


respectively  (see  Crimer  (1964)).  Here  vk  is  the  kth  moment  of  the 

distribution  in  question.  Suppose  we  look  at  the  coefficient  of  1/n 

-x 

in  (  8)  for  the  case  of  an  exponential  distribution  with  density  e 

for  x  >  0.  It  is  195  and  in  terms  of  this  the  implication  is  that  one 

would  need  a  sample  size  of  about  600  to  get  a  variance  of  the  order 

of  magnitude  of  one.  The  case  is  much  more  extreme  for  the  coefficient 

^  |  x  |  1 

of  1/n  in  (9)  for  the  case  of  a  symmetric  exponential  density  e  y 
The  coefficient  is  39,744. 

Deconvolution  weights.  Here  we  will  sketch  an  argument  that  allows  us  to 
get  an  asymptotic  approximation  for  the  covariances  of  the  principal 
random  part  of  deconvolution  weight  estimates  b^.  A  similar  argument 
can  be  used  to  obtain  such  an  approximation  for  the  covariances  of  the 
principal  random  part  of  the  estimates  ak.  Expression  (  7 )  can  be  rewritten  as 
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J/2 

X 


_JL  G  00 

b.  -  TirZ-/  (2*  fa,)]  1  COS  { -G  (X . )  +  X.  +  kX  ) 

k  J+l  n  j  n  j  *  j  j 


Now 


.  x  I  f  (xj-k  a4J\" 

(10)  (f^Xj))'  *-  (E  f.a,))-*  +  nEj'<x')  J-j 

i  /  ,  f  (X  )-Ef  (X  ) 

■  <E  (X  -  T  " *  +  °(fn<V 


Efn<V> 


Further 


(11)  cos 


G  (it) 

-Gn(Xj>  +“V-XJ  +kXJ 


EG  (it) 


-  cos  {-  EG  (X x)  +  — - - X,  +  kX. 

n  j  *  j  j 


+  [-  W  +  *W  +- 


G  Or) -EG  Or) 


n 


J| 


EG  Or) 

-  cos  {-  EG  (X.)  +  - - X.  +  kX. 

n  j  ”  J  J 


-  sin 


-  EG  (X.)  + 
n  J 


EG  (ir ) 


- - X.  +  kX.. 

*  J  J] 


r  G  (ir)-EG  (ir) 

[-  Gn(XJ>  +  *.<XJ)  + - 1 - X 


n  ■  .  J 


r  G  (ff)-EG  (.)  1 

+  °pL‘  Gn<xj> +  EG„<V  +  — ; - XjJ  • 


First  it  should  be  noted  thst  the  second  tern  on  the  right  of  (10)  will  be 
of  smeller  order  then  the  second  term  on  the  right  of  (LI).  This  implies 
that  the  prlncipel  random  pert  of  b^  (the  deterministic  mean  is  neglected 
here)  can  be  approximated  by 


(2ir  ftXj))”1  sin  (-  h^Xj)  +  kXj) 
r  G  (ir)-EG  (ir)  1 

[-  W  +  BGn(V  "  XjJ 


The  principal  part  of 


G  00 -EG  (ir) 

-  G  (X)  +  EG  (X)  ♦  -2 - - - X,  0  <  X  <  v  , 

n  n  n 


asymptotically  has  the  covariance  (the  argument  is  like  that  given  for 
Rn(x)  in  (  4)) 


4tt 


a4  2 
A  ny 


|min  y  W2(u,v,w)  dudvdw 


if  A  n  -*•  «,  A(n)  0.  This  implies  that  the  covariance  of  principal 

random  parts  of  ,b^  are  (j ,k  fixed) 


f'f 

‘2  JJ 


An y  *'0 


(f(X)f(p))”*  sin  (h^X)  -  kX) 


sin  (h^vO-jp)  jmin  2-j  -  j  dXdy 

y ^2cu.v,i 


,w)  dudvdw  . 
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